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ABSTRACT 

The metallicities implied by collisionally excited lines (CELs) of heavy elements in H II regions 
are systematically lower than those implied by recombination lines (RLs) by factors ~ 2, introducing 
uncertainties of the same order in the metallicities inferred for the interstellar medium of any star- 
forming galaxy. Most explanations of this discrepancy are based on the different sensitivities of 
CELs and RLs to electron temperature, and invoke either some extra heating mechanism producing 
temperature fluctuations in the ionized region or the addition of cold gas in metal-rich inclusions or 
ionized by cosmic rays or X rays. These explanations will change the temperature structure of the 
ionized gas from the one predicted by simple photoionization models and, depending on which one is 
correct, will imply different metallicities for the emitting gas. We select nine H II regions with observed 
spectra of high quality and show that simple models with metallicities close to the ones implied by 
oxygen CELs reproduce easily their temperature structure, measured with T e ([N II])/T c ([0 III]), and 
their oxygen CELs emission. We discuss the strong constraints that this agreement places on the 
possible explanations of the discrepancy and suggest that the simplest explanation, namely errors in 
the line recombination coefficients by factors ~ 2, might be the correct one. In such case, CELs will 
provide the best estimates of metallicity. 
Subject headings: ISM: abundances — H II regions 



1. INTRODUCTION 

In an H II region, the temperature of a volume of 
gas in thermal equilibrium is determined by the balance 
between heating, mainly through photoionization, and 
cooling, mainly through recombination and radiation of 
collisionally excited lines (CELs). The most important 
coolants are CELs arising from low-energy levels of ions 
of abundant metals, like + , ++ , N + , S + , and S ++ , 
and the gas metallicity is the single most important pa- 
rameter in the determination of the final temperature of 
the gas (|Osterbrock fe Ferlandl[2r3ol ). 

Temperatures can be measured using the relative in- 
tensities of CELs emitted by ions like ++ and N + , and 
their associated temperatures, T c ([0 III]) and T e ([N II]), 
can be used to characterize the zones of high and low 
degree of ionization within a nebula. These temper- 
atures can be used to calculate line emissivities and 
hence the abundances of different elements. The oxy- 
gen abundance is probably the most reliable one that 
can be derived in H II regions and it is generally used 
as a proxy for the metallicity of the gas. When optical 
spectra of high quality are available, the oxygen abun- 
dance is derived by adding the + /H + and ++ /H + 
abundance ratios implied by the intensities of optical 
[O II] and [O III] lines (relative to hydrogen recombi- 
nation lines) and the corresponding emissivities calcu- 



lated for T ([N II]) and T c ([0 III]), respectively. How- 
ever, the ++ abundance can also be determined us- 
ing several weak optical O II recombination lines (RLs) 
and these lead to abundances that are consistently higher 
than the ones derived from CELs. A similar result has 
been found with CELs and RLs of other io ns, both in H II 
regio n s and in planetar y nebulae (see, e.g. JEsteban et "all 
l200i iLiu et all l200l . but with higher uncertainties. 
The abundance discrepancy factors (ADFs) for ++ , 
ADF{0 ++ ) = (0 ++ /H+) rl /(0 ++ /H+)cel, are of the 
order of 2 for H II regions and many planetary nebulae, 
but can reach muc h higher values, up to 70, in the latter 
objects (|Liul 12006. and references therein). 

Most of the proposed explanations of the ADF are 
based on the different dependence with temperature of 
the emissivities of CELs (strongly increasing with T e ) 
and RLs (oc T~ m with m ~ 1). If there are large temper- 
ature fluctuations in the gas, the emission of CELs will 
be heavily weighted towards the hotter regions; the emis- 
sion of RLs towards the colder ones. Depending on the 
mechanism responsible for the fluctuations, either CELs 
or RLs will provide the better estimate of the gas metal- 
licity. The fluctuations could be due to hot gas heated 
by shocks, stell ar winds, small dust grains, magnet ic re- 
connection (see iTorres-Peimbert fc Peimb ert 20031 an d 
references t herein) or can be due to cold gas in metal-rich 



Electronic address: mrodri@inaoep.mx 
Electronic address: jogarciaSSiac.es 

1 Visiting Astronomer, Instituto Nacional de Astrofi'sica, Optica 
y Electronica, Apdo Postal 51 y 216, 72000 Puebla, Mexico 

2 Current address: Instituto de Astrofi'sica de Canarias, E-38200 
La Laguna, Tenerife, Spain 



inclusions (ILiu et all l20"o"(it iTsamis fc PeauignotJ 120051 : 
Stasihsk a et al.l 2007 ) or cold gas ionized b 



b y cosmic 

rays (|Giammanco fc Bec kman 2 0051) or X rays (|Ercolanol 
2009). Some or all of these effects are likely to be present 
in real objects, but it is not clear yet whether they can 
explain the derived ADFs or how much they contribute 
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to the observed emission. This introduces an uncertainty 
in the absolute values of the metallicities in the interstel- 
lar medium. 

One way to approach this issue, the one we will 
follow here, is to build a series of simple photoionization 
models with input parameters that can be considered 
representative of some observed objects and to see 
whether the models reproduce the observed spectra. We 
center on nine H II regio ns with available sp ectra of high 
quality and use Cloudy (jFerland et al.l[l99a ) to calculate 
grids of photoionization models that depend only on a 
few free parameters. The observed and predicted spectra 
are analyzed in a similar way to determine electron 
densities, temperatures, and the ionic and total abun- 
dances of oxygen using the intensities of CELs relative 
to hydrogen lines. We show that the models reproduce 
easily the general trends defined by the observations. 
We argue that they are representative of the observed 
objects and discuss how far they can be modified by 
adding some extra ingredient that could also explain the 
intensities of O II RLs without destroying the agreement. 

2. THE OBSERVATIONAL SAMPLE 

We consider nine H II regions with published spec- 
tra of the best quali ty: the Galacti c H II regions M42 
(lEsteban et al.l 12004ft. NGC 3576 (iGarcia-Roias et all 
12004ft. NGC~ 2467 flOarria-Roias et al. || 2005ft . M16. M20. 
NGC 3603 (IGarcia-Roias et all l2006ft~~~M8 and M17 
(IGarcia-Roias et al.l 12007ft . and 30 Doradus (|Peimbertl 
12003ft in the Large Magellanic Cloud (LMC). Deep spec- 
tra of these nebulae were obtained with the VLT UVES 
echelle spectrograph covering the range 3100-10400 A 
with spectral resolution of A A ~ 8800/ A. The spectra, 
which were extracted over small fractional areas, 3" x 8'.'5 
or 3" x 10", allowed the measurement of 235-555 emis- 
sion lines in each nebula, including several weak O II 
ORLs, that were detected with good signal-to-noise ra- 
tios. Physical conditions could be derived from several 
diagnostics that led to broadly consistent values. 

Since we will be comparing the values of the tempera- 
tures and abundances derived from the observed spectra 
with those derived from the line intensities predicted by 
the photoionization models, we decided to use the same 
atomic data for all calculations. Hence, we recomputed 
the parameters of interest with the nebular package in 
IRAF 3 but using the same atomic data as Cloudy. This 
implied changing the tables of transition probabilities 
and collision strengths used in IRAF for the elements 
of interest (references for th e values we used for thes e 
atomic data can be found in IGarcia-Roias et alj [2009ft . 
We did not change the number of levels used in the cal- 
culations nor the procedure used by IRAF to get the 
emissivities of the hydrogen lines, but the effects of these 
differences are small. 

The recomputed electron densities, tIq , are a 
weighted mean of the values implied by the diagnostics 
[S II] A6716/A6731, [O II] A3726/A3729, and [CI III] 
A5517/A5537. The values of T e ([0 III]) are based on the 
ratio of line intensities (A4959 + A5007)/A4363; those of 
T e ([N II]) on (A6548 + A6583)/A5755. In some of the 

3 IRAF is distributed by NOAO, which is operated by AURA, 
Inc., under cooperative agreement with NSF. 



previous studies, the values of T e ([N II]) were corrected 
for the contri bution of rec ombination to the intensity of 
[N II] A5755 (|Rubinl 1 1986ft . but since the corrections are 
based on uncertain estimates of the abundance of N ++ , 
we prefer to apply that correction to the models (see 
Section [3J>. The values of J([0 II] A3726+29)//(H/3), 
n e , and T e ([N II]) are used to derive the + /H + abun- 
dance ratio; /([O III] A4959 + A5007)/J(H/3), n e , and 
T e ([0 III]) are used for ++ /H + . The total abundance 
of oxygen is then obtained from the sum of the + and 
++ ionic abundances. We also recomputed the ++ 
abundances implied by RLs using the O II multiplet Ml 
and t he recombination coefficients of IPeauignot et al.l 
(|1991ft which are those implemented in Cloudy. In 
Table [1] we show the physical conditions, T e , n e and the 
degree of ionization, the O/H abundances derived from 
CELs and the ADFs for all the sample objects. Despite 
the slight changes in the procedure and in the atomic 
data used, the final values of the oxygen abundance and 
of the ADF are consistent within the errors with the 
values derived in the previous analyses. 

3. THE PHOTOIONIZATION MODELS 

Calculations were perfo rmed with ve r sion 8.00 of 
Cloudy, last described by iFerland et ail (| 1998ft . Tak- 
ing into account that the observations covered small 
fractional areas of the objects, that different ions lead 
to similar densities, and that the line emissivities are 
roughly proportional to the square of the density, one- 
dimensional plane-parallel models with constant density 
should be a good approximation. 

We use models with total hydrogen densities close 
to the electron densities derived for the observed ob- 
jects: n H = 300 cm" 3 for M17, M20, 30 Doradus, and 
NGC 2467; n H = 1000 cm" 3 for NGC 3576, M8, and 
M16; and n H = 5000 cm-' 3 for NGC 3603 and M42. 

We define the ionizing radiation field through two pa- 
rameters: the effective temperature of the ionizing star, 
T e ff , and the ionization parameter u — 4>n/nnc, with <fin 
the number of hydrogen ionizing photons arriving at the 
inner face of the cloud per unit area and unit time, and 
c the speed of light. 

We consider the two sets of state-of-the-art model 
atmospheres for O-type sta rs incl uded in Cloudy: 
WM-BASIC (IPauldrach et all 120011 ) and TLUSTY 
(jLanz &; Hubenvl 12003ft . Comparisons between the 
ionizing fluxes predicted by different stellar models and 
discussions on the effects on t he resulting photoion - 
ization models can be fo und i n iMartins et al.l (2005), 
iSimon-Diaz fc Stasihska (2008), or in the references 
therein. Here we present the results obtained with 
wm-basic and use the results obtained with TLUSTY 
as a check for robustness. We use stellar models with 
surface gravity log(g) = 4.0, solar metallicity, and 
effective temperatures T cS = 35000, 37000, 39000, 
41000, 43000, 45000, and 50000 K. The stellar models 
with T cff = 35000 K to 45000 K represent roughly th e 
spectral types from 03 to 08 (jMartins et al.l [2005); 
the models with T c g = 50000 K probe the effects of 
more massive stars. The values of u are chosen for 
each value of T e g so that the set of models cover the 
range of degrees of ionization found for the observed 
objects. For the models with nn = 300 and 1000 
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TABLE 1 

Physical conditions and oxygen abundances for the sample objects 





n c 


T„(TN III) 

-teU 1 ^ - LJ -J/ 


tjio hid 












Object 


(cm -3 ) 


(K) 


(K) 


12 + log 


(O/H) 


log(0+/0++) 


ADF{0++) 


M8 


1550 ± 150 


8500 ± 120 


8020 ± 100 


8.44 


± 


0.03 


+0.40 ± 0.05 


2.2 + 0.2 


M16 


980 ± 120 


8500 ± 150 


7580 ± 150 


8.52 


± 


0.03 


+0.45 ± 0.06 


2.4 + 0.6 


M17 


420 ± 80 


9150 ± 250 


7950 ± 100 


8.53 


± 


0.02 


-0.76 ± 0.06 


1.8 + 0.2 


M20 


220 ± 40 


8500 ± 150 


7750 ± 250 


8.48 


± 


0.03 


+0.64 + 0.07 


1.9 ± 1.2 


M42 


6800 ± 600 


10100 ± 250 


8250 ± 60 


8.53 


± 


0.02 


-0.68 + 0.05 


1.50 + 0.05 


NGC 2467 


260 ± 50 


9600 ± 200 


8900 ± 100 


8.35 


± 


0.03 


+0.35 + 0.05 


1.8 + 0.2 


NGC 3576 


1400 ± 200 


8950 ± 200 


8400 ± 80 


8.53 


± 


0.02 


-0.40 ± 0.06 


2.0 + 0.3 


NGC 3603 


2350 ± 400 


11650 ±570 


9000 ± 150 


8.47 


± 


0.03 


-1.21 + 0.09 


1.9 + 0.3 


30 Doradus 


380 ± 50 


10800 ± 300 


9850 ± 100 


8.36 


± 


0.02 


-0.77 + 0.06 


1.5 + 0.1 



cm 3 we use log(u) = —1.5, —1.0, —0.5, and 0.0 for 
T cS = 35000 K; log(u) = -2.5, -2.0, -1.5, and -1.0 for 
T cff = 37000 K; log(u) = -3.0, -2.5, -2.0, and -1.5 for 
T cff = 39000 K; and log(u) = -3.5, -3.0, -2.5 and -2.0 
for T e $ = 41000, 43000, 45000, and 50000 K. For the 
models with nn = 5000 cm -3 , we increased the values 
of log(u) by 0.5 dex. 

We also computed models ionized by two stars with 
T c s = 35000 and 45000 K that provide equal amounts of 
hydrogen ionizing photons (log(u) = —3.3, —2.8, —2.3, 
or —1.8 for each star). These models illustrate the effect 
of a composite stellar radiation field. 

We use the set of abundances labeled as "H II region" 
(or "Orion" ) in Cloudy and the associated "H II region" 
dust grains. We scaled both metals and grains in the 
"H II region" set by the same factor in order to explore 
different metallicities. We changed the metallicity using 
steps of 0.1 dex and explored the range 12 + log(0/H) = 
8.3-9.1. 

By default, Cloudy stops when the model reaches an 
electron temperature of 4000 K. However, since some of 
the models with high metallicity reach this temperature 
in the fully ionized zone, we switched off this criterion 
and used instead the condition log(n e /nn) < —0.5. 

For each model we get a list with the intensities relative 
to H(3 of the lines of interest. We introduce this list in 
IRAF and use it to calculate physical conditions and ionic 
abundances following the same procedure we used in the 
analysis of the observed objects (see Section [2]). 

As mentioned in Section [2j the upper level of the 
[N II] A5755 line used to derive TJ\ N II}) can b e 
populated by recombination processes (Rubin 1986). 
If this effect is not taken into account, the values of 
T e ([N II]) derived for the models will be lower than those 
derived for the observed objects, especi ally for th e cases 
where the degree of ionization is high. iLiu et aLl (2000) 
estimated the contribution to the intensity of this line in 
terms of the N ++ abundance. We used their correction 
to add this contribution to the intensity of the [N II] 
in each model using the N ++ and H + column densities 
and T c = 10000 K (the dependence on T is very weak). 
Then we calculated a second set of physical conditions 
and abundances for all the models. Since the correction 
is just an estimate, both the corrected and uncorrected 
results are shown below. 

4. RESULTS 

We compared the physical conditions and abundances 
derived for the models with those derived for the ob- 



served objects. Most of the models lead to derived abun- 
dances, (0/H) out , that are l ower than the i nput ones. 
This is a well known effect (Stasihska 2005 and refer- 
ences therein), due to the exponential dependence with 
temperature of the CELs used to derive T e and the + 
and ++ abundances. In the presence of a temperature 
gradient, this dependence leads to an overestimation of 
T e and hence of the emissivity of the CELs. This trans- 
lates into an underestimation of the d erived abundance s 
relative to hydrogen. As discussed by iStasihskal (|2005l) . 
under some circumstances, models with high metallic- 
ity and low density develop a gradient so steep that 
the difference between the input (O/H) and (0/H) out 
becomes substantial. One extreme case is our model 
with n H = 300 cm" 3 , T off = 41000 K, log(u) = -2.0, 
and 12 + log(0/H) = 9.1. The abundance derived from 
the analysis of the spectrum predicted for this model is 
12 + log(0/H) out = 8.55, and ADF(0 ++ ) ~ 4. Other 
high-metallicity models have ADFs close to the ones de- 
rived for the observed objects, but also T c ([0 III]) < 6000 
K, much lower than the observed values. All the models 
that consistently give values of (0/H) out , T C ([N II]), and 
T c ([0 III]) close to the observed ones, have input metal- 
licities to 0.06 dex higher than the derived ones and 
ADF(0 ++ ) ~ 1.1. 

Figures [l] and [2] show a comparison of T c ( [O III] ) and 
T e ([N irj)/T c ([0 III]) as a function of the derived val- 
ues for the degree of ionization, log(0 + /0 ++ ), for mod- 
els and observations. The results implied by models 
with two input metallicities are shown for each object: 
12 + log(0/H) = 8.5 and 8.6 for M8, M16, M17, M20, 
NGC 3576, M42, and NGC 3603; 12 + log(0/H) = 8.3 
and 8.4 for NGC 2467 and 30 Doradus. The output 
metallicities of these models, (0/H) out , bracket the ob- 
served values listed in Table [TJ The observed objects are 
compared with those models that have similar values of 
n c . It can be seen that these models also reproduce eas- 
ily the temperatures derived for the observed nebulae. 
An exception is the value of T C ([N II])/T e ([0 III]) for 
NGC 3576. This could be just a statistical fluctuation, 
with one in nine measurements deviating with a 2a er- 
ror bar from the expected result, or could indicate that 
something is different in the environment of this object. 

Note that the models ionized by two stars with T c g = 
35000 K and 45000 K (the pentagons in Figures Q] and [2]) 
give results that are intermediate between those shown 
by the models ionized by one star with either T ff = 
35000 K or 45000 K (the asterisks and the circles, re- 
spectively). Hence, models ionized by composite spectra, 
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Fig. 1. — Values of T c ([0 III]) and T e ([N II])/T c ([0 III]) for observed objects (filled black squares) and photoionization models 
(other symbols). The models have rijj = 300 cm -3 and T c g = 35000 K (red asterisks), 37000 K (orange open squares), 39000 K (green 
six-pointed stars), 41000 K (cyan triangles), 43000 K (blue four-pointed stars), 45000 K (violet circles), 50000 K (black inverted Ys), or 
both T c fj = 35000 K and 45000 K (magenta pentagons). Big symbols correspond to models with 12 + log(0/H) = 8.6 (top panels) and 
8.4 (bottom panels). Small symbols have 12 + log(0/H) = 8.5 (top panels) and 8.3 (bottom panels). Dotted symbols show values derived 
considering the contribution of recombination to the intensity of [N II] A5755. The small filled circles show the effects of intervening 
material changing the shape of the ionizing radiation field (top panels, for models with T c ff = 39000 K and 12 + log(0/H) = 8.6) and of 
decreasing the abundances of C, N, S, and Ne relative to O (bottom panels, for models with T c g = 45000 K and 12 + log(0/H) = 8.4). 
See the text for more explanations. 



with several values of T e g, will also appear in the regions 
sampled by the single T c s models in Figures Q] and [2] 

Very similar results to the ones found with WM-BASIC 
were obtained using the tlusty model atmospheres. Al- 
though individual models would appear in somewhat dif- 
ferent positions if plotted in Figures [T] and [2l the areas 
covered by both kinds of models are roughly similar. 

Many simplifying assumptions were made constructing 
these models. The agreement with the observations sug- 
gests that despite the simplifications and uncertainties, 
the models have captured some basic characteristics of 
the observed objects. We have already addressed the ef- 
fects of having a composite ionizing spectra and of using 
different model atmospheres. We now address other un- 
certainties and simplifications of the models to further 
explore how general are the results. 

4.1. Effect of intervening material 

Our models assume that there is no intervening ma- 
terial between the ionizing stars and the observed area. 
However, this material, if present, will change the shape 
(and luminosity) of the ionizing radiation field reaching 
the observed region and, in principle, this could lead to 



departures from the behavior of the models in Figures [T] 
and O To check for this effect we calculated several 
models where the stellar radiation field was first passed 
through a slab of material of fixed depth and the trans- 
mitted radiation was then used to ionize the slab of in- 
terest. 

We start considering intervening material with nn = 
300 cm" 3 , 12 + log(0/H) = 8.6, and a radiation field 
with T cff = 39000 K and log(u) = -1.0. If enough mate- 
rial is present, these conditions will produce a depth of 
ionized material of about 4 pc. We consider slabs with 
thicknesses of 1, 2, and 3 pc and use the radiation trans- 
mitted through them to ionize material with the same 
riH and metallicity. These models have values of the ion- 
ization parameter of log(u) = —1.23, —1.5, and —1.93, 
respectively. The spectra produced by the latter slabs 
were analyzed in the same manner as the original mod- 
els and the results are plotted as small filled circles in 
the top panels of Figure [TJ It can be seen that they give 
results similar to the original grid of models. 

We calculated similar models (with 12 + log(0/H) = 
8.6, T cff = 39000 K, and initial log(u) = -1.0) for ma- 
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Fig. 2. — Values of T e ([0 III]) and T C ([N II])/T c ([0 III]) for observed objects (filled black squares) and photoionization models (other 
symbols). The models have njj = 1000 cm -3 (top panels) or 5000 cm -3 (bottom panels), and T c ff = 35000 K (red asterisks), 37000 K 
(orange open squares), 39000 K (green six-pointed stars), 41000 K (cyan triangles), 43000 K (blue four-pointed stars), 45000 K (violet 
circles), 50000 K (black inverted Ys), or both T e g = 35000 K and 45000 K (magenta pentagons). Big symbols correspond to models 
with 12 + log(0/H) = 8.6; small symbols have 12 + log(0/H) = 8.5. Dotted symbols show values derived considering the contribution 
of recombination to the intensity of [N II] A5755. The small filled circles show the effects of decreasing the metallicity of the ionizing 
stars by a factor of 2 (top panels) and of changing the dust grains from "Orion" type to "ISM" type (bottom panels) for models with 
T cff = 39000 K and 12 + log(0/H) = 8.6. Models with T eff = 39000 K and 12 + log(0/H) = 8.6 but with no grains included are shown as 
crosses in the bottom panels. See the text for more explanations. 



terial with njj = 1000 cm -3 . The incident radiation 
field was first passed through slabs with the same den- 
sity and depths of 0.3, 0.6, and 0.9 pc (so that for the 
final models log(w) = -1.23, -1.5, and -1.9). We did 
other tests with intervening material of lower density, 
7iH = 100 cm~ 3 , and depths of 1, 2, and 3 pc (leading 
to log(u) = —1.06, —1.11, and —1.17 for the final mod- 
els with ripr = 1000 cm" 3 ). All these models also gave 
similar results to those produced by the original grid of 
models. 

4.2. Changing the relative abundances of some elements 

We are using the set of "H II region" chemical 
abundances in Cloudy (based on the abundances de- 
rived for M42 bv iBaldw in et al.|[l99l iRubin et~aT1ll991l; 
IQsterbrock et al.lll992T ) and scaling all of them by the 
same amount to change the overall metallicity of the 
gas, which we measure through O/H. For the other el- 
ements, not all ionization states show lines in the op- 
tical and their total abundances are derived using ion- 
ization correction factors that make them less reliable. 
Changing the abundances relative to oxygen of those el- 
ements that contribute significantly to the cooling might 



change our results. We compared the values of the C/O, 
N/O, S/O, and Ne/O abundance ratios in Cloudy with 
the values derived for these ratios in the original papers 
where the analysis of the observed objects was first pre- 
sented. The most discrepant results are found for 30 Do- 
radus, with C/O, N/O, S/O, and Ne/O lower than the 
Cloudy values by 0.4, 0.5, 0.1, and 0.14 dex, respectively. 
Hence, we ran three new models with njj = 300 cm~ 3 , 
12+log(0/H) = 8.4, T cff = 45000 K, log(u) = -3.0 -2.5, 
and —2.0, and with the abundances of C, N, S, and Ne 
scaled down by the amounts listed above. The results 
are shown as small filled circles in the bottom panels of 
Figure [TJ As expected, they show somewhat higher val- 
ues of T e , but the values of T e ([N II])/T ([0 III]) do not 
change significantly. 

4.3. Changing the metallicity of the ionizing star 

All the calculations were performed with stellar spec- 



tra f or solar metallici ty, with 12 + log(0/H) f 



> 



.7 



(e.g. iScott et all 120091 and references therein). How- 
ever, the models have input metallicitics in the range 
12 + log(O/H) = 8.3-8.6. We calculated four new mod- 
els with n H = 1000 cm" 3 and 12 + log(0/H) = 8.6, ion- 
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ized by stars with T c g = 39000 K and half-solar metallic- 
ity, and with log(u) = -3.0 -2.5, -2.0, and -1.5. The 
results of these models are plotted as small filled circles 
in the top panels of Figure [2] The new models have a 
higher degree of ionization but show a similar behavior 
to the original ones. 

4.4. Effect of dust grains 

The models include dust grains, and the grains affect 
the temperature, mainly through heating by the photo- 
electric effect. The fractional contribution of grains to 
the heating increases with the ionization parameter and 
decreases with the effective temperature of the ionizing 
star. For the main photoionization models described so 
far, the contribution can be negligible or reach ~ 30% 
of the global heating. We are using the "H II region" 
(or "Orion") dust grains, a mixture of graphite and sil- 
icate grains that is deficient in small pa rticles in order 
to ex plain the grey extinction in M42 (Baldwin et al. 
1991), but there are many uncertainties related to which 
grains would be suitable. Changing the type of dust 
grains will change the temperature structure. To illus- 
trate this effect, we calculated four new models with 
n H = 5000 cm" 3 , 12 + log(0/H) = 8.6, T cS = 39000 K, 
and log(u) = —2.5, —2.0, —1.5, and —1.0, changing the 
type of dust from "H II region" to "ISM" (interstellar 
medium). The "ISM" dust grains are also composed of 
graphite and silicate, but have a size distribution that 
includes small grains. The results of these models are 
plotted as small filled circles in the bottom panels of 
Figure [3 They have higher values of T c ([0 III]) and 
lower values of T e ([N II])/T e ([0 III]). Whereas the origi- 
nal models (with T cS = 39000 K, 12 + log(0/H) = 8.6, 
and "H II region" dust) had 8% as the maximum contri- 
bution of grains to the global heating, the new models 
reach 30%. Note however, that the results are not very 
different from the original ones. 

The same test models were computed with no grains 
included (but keeping the depleted abundances of re- 
fractory elements). These models show lower values for 
the derived electron temperatures and higher values of 
T C ([N II])/T e ([0 III]) than the original ones. This is 
shown by the crosses in the bottom panels of Figure [2J 

One of our objects, 30 Doradus, is located in the LMC, 
where dust grains can be very different from Galactic 
grains. However, the amount of dust in 30 Doradus is 
probably well reproduced by the comparison models of 
the bottom panels in Figure Q] All our models have high 
depletion factors of refractory elements onto dust grains, 
and the depletion factors are kept constant by scaling 
together the abundances of metals and grains. If we con- 
sider iron a good representative of refractory elements, 
this assumption is supported by the high iron depletion 
factors found in 30 Doradus and several Galactic H II 
regions ((Rodriguez fc Rubinl 120051 ). On the other hand, 
the grains in 30 Doradus could have a different composi- 
tion and be smaller than their Galact ic counterparts (see 
ISloan et al.|[200l IParadis et al.ll2009L for some recent re- 
sults), but the effects of the different characteristics of 
grains in 30 Doradus on the calculated temperatures are 
likely to be of the order of those discussed above. For a 
more thorough investigation of the effects of changes in 
composition and size distribution on model output see 
Ivan Hoof etall (|200l . 



4.5. Advection 

The models we are calculating are static. In the 
dynamical case, the ionization of neutral material 
entering the ionization front (advection), can in- 
crease the temperature of the gas in th at region (e.g. 
iRodriguez-Gaspar fc Tenorio-Tagld 1 19981) . This will 
translate into higher values of T fi ([N Hj), but the effe ct 
is probably small, as discussed bv lHennev et a l. (2005). 



5. RESULTS FOR THE CASE WITH EXTRA HEATING 

We have shown that the grids of simple photoioniza- 
tion models reproduce very well the temperatures and 
metallicities implied by the intensities of CELs relative 
to hydrogen recombination lines in the observed objects. 
However, these models have ADF(0 ++ ) ~ 1.1 and hence 
do not reproduce the intensities of O II RLs. In order to 
reproduce the intensities of both CELs and RLs, some 
extra ingredient must be added to the models. What 
characteristics should this extra ingredient have so that 
it does not destroy the agreement found so far between 
observations and models? 

We start exploring the case where the extra ingredient 
is a heating agent that acts in some regions of the neb- 
ula producing temperature fluctuations. In this case, the 
emission of CELs will be strongly biased towards the hot- 
ter regions leading to overestimates of the derived tem- 
peratures and underestimates of the abundances. The 
abundances derived using O II RLs and the tempera- 
ture f luctuations formalism first develo ped by iPeimbertl 
(|1967f ) (see also iPeimbert et al.l 12004 and references 
therein) will be more reliable. As an example, in M8 
and M16 the oxygen abundances de rived with this pro- 
cedure are 12 + log(0/H) ~ 8.7-8.8 (|Garcia-Roias et al.l 
20061 12001 . 

We modified the function that can be used to introduce 
a term of extra heating in Cloudy so that the term was 
only added at localized depths within the ionized slab. 
In this way we can produce spikes in the temperature 
structure of the ionized region: the spikes appear in those 
zones where the term of extra heating is different from 
zer o. The approach is s omewhat similar to the one used 
by iBinette et al.l ((2001). The temperature fluctuations 
due to these spikes introduce a bias in the temperatures 
and abundances derived from CELs that produce in turn 
values of ADF(0 ++ ) significantly larger than 1. Note 
that as long as the spikes are regularly distributed across 
the ionized region, models with many thin spikes lead 
to results similar to those obtained from models where 
these spikes are grouped producing fewer, but broader, 
spikes. If one fixes the total number of spikes appear- 
ing in the ionized slab, the required values for the ADF 
can be obtained by varying the amount of extra heating 
through changes in the width or the height of the spikes. 

The upper panel of Figure [3] shows the temperature 
structure of two models that produce ADF(0 ++ ) ~ 2 for 
some of the conditions suitable to describe M8 or M16. 
The middle and lower panels of Figure [3] show for these 
two models the emissivities of [O III] A5007 and multiplet 
1 of O II with arbitrary normalization. Since the models 
have almost constant electron densities, the emissivities 
are modulated by the ionization fraction of ++ and by 
the electron temperature. These models illustrate how 
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Fig. 3. — Upper panel: Values of T e for two models with terms of 
extra heating added in selected regions, as a function of normalized 
depth into the ionized slab (with depths of 2.49 X 10 17 cm for the 
model plotted with a solid line and 2.66 X 10 17 cm for the model 
plotted with a dotted line). Both models have T c g = 39000 K, 
log(w) = -2.5, n H = 1000 cm" 3 , and input 12 + log(0/H) = 8.75. 
Both models produce ADF(0 ++ ) ~ 2. Middle and lower panels: 
Emissivities of [O III] A5007 and multiplct 1 of O II for the same 
models in arbitrary units. 



large the temperature fluctuations must be in order to 
produce the required ADFs. The models have T e s = 
39000 K, log(u) = -2.5, n H = 1000 cm" 3 , and input 12+ 
log(0/H) = 8.75. Derived values are: log(0+/0 ++ ) ~ 
0.4, 12 + log(0/H) = 8.56 and 8.60 (for the models with 



the broad and thin spikes, respectively), and T c ( 



9300 K, T c ([0 III 
8000 K, T c ([0 III 



NII]) = 
NII]) = 



9800 K (broad spikes), T e ( 
8200 K (thin spikes). In the model 
with broad spikes, the extra heating contributes 32% of 
the model's total heating; in the model with thin spikes 
this contribution is just 6%. 

Models with a term of extra heating just 20% lower 
than the one used for the model with broad spikes in 
Figure [3] do not have enough temperature fluctuations 
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Fig. 4.— Values of T e ([N II])/T e ([0 III]) as a function of 
ADF(0++) for models with T cff = 39000 K, log(u) = -2.5, 
Tin = 1000 cm -3 , input 12 + log(0/H) = 8.75, and temperature 
spikes of different thickness and widths (four-pointed stars). Also 
plotted are those models from the top panels in Figure [2] that 
have log(0+/0"'"+) ~ 0.45 (open circles), a value similar to those 
shown by M8 and M16 (plotted with filled squares). 



to reach ADF(0 ++ ) - 2: they have ADF(0 ++ ) < 1.8 
for spikes of any width. On the other hand, we did not 
explore models with terms of extra heating larger than 
the one used for the model with thin spikes in Figure [3j 
Those models would require very short integration steps 
in Cloudy. Besides, the thin spikes would dissipate very 
quickly, as discussed below. 

The result found for the two models in Figure [31 
T e ([N II]) < T c ([0 III]), is characteristic of all the 
models that can be constructed in this way to produce 
large ADFs. This can be seen in Figure [3J where we 
plot the values of T ([N II])/T ([0 III]) as a function of 
ADF(0 ++ ) for several models with temperature spikes 
of different thicknesses and widths. Also shown are 
the values derived for the original models with no ex- 
tra heating and plotted in the top panels of Figure [2] 
with log(0 + /0 ++ ) ~ 0.45 (a similar degree of ionization 
to the ones found in M8 and M16). None of the mod- 
els with extra heating can reproduce the observed values 
of the temperature ratio, both because the extra heat- 
ing tends to make similar the average temperatures of 
the [N II] and [O III] emitting regions, and because the 
[O III] diagnostic ratio is more sensitive to increments 
in temperature than the [N II] ratio. In order to get 
T C ([N II]) > T e ([0 III]) as in the observed objects, the 
term of extra heating would have to increase with depth 
into the cloud. This explanation of the ADF would also 
work if the temperature spikes are created continuously 
across the full ionized region but dissipate more easily 
in the regions closer to the star. Any attempt to ex- 
plain the abundance discrepancy with temperature fluc- 
tuations due to extra heating would need to explain this 
effect in a natural way. 

The results described in this section are based on 
illustrative but unrealistic models. The pressure differ- 
ences introduced by the high temperature regions would 
produce instabilities, and pressure equilibrium would 
be reached in times given by the ratio of the width 
of the spikes to the speed of sound in the gas (~ 10 
km s" 1 ). Hot regions with widths of 10 15 cm, like the 
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ones in the model with thin spikes in Figure [31 would 
have dynamical lifetimes of ~ 30 years. Heating and 
cooling processes have simila r timescales in ionized gas 
(jOsterbrock fc Ferlandl 120061 1 . Hence, the mechanism 
responsible for the hot spikes should renew th em within 
this timescale, as previously noted by I Ferlandl (|200I[ ). 

6. COMMENTS ON MODELS THAT INCLUDE COLD GAS 

Models that reproduce the intensities of O II RLs 
can also be built by adding a component of cold gas 
that contributes significantly to the emission of RLs but 
not to the emission of CELs. This cold gas could arise 
from the photoionization of metal-rich inclusi ons (see 
iTsamis fc Pequignotll2005HStasihska et alj|2007l for H II 
regio ns) or could be g as ionized by X ravs (lErcolanol 
2009) or cosmic rays (jGiammanco fc Beckmanl 120051) . 
Because the cold gas is poor in hydrogen or because 
its ionization fraction of ++ is larger than the one for 
H + , CELs will give better estimates of the metallicity in 
the models with cold gas that have consid ered this issue 
(see iStasihska et~aT1 l2007t lErcolanol 12009( 1 . To explain 
the ADF, the gas should have a significant fraction of O 
as ++ and a temperature low enough to suppress the 
emission of CELs. Some of the temperatures considered 
in the works mentioned above for the cold component 
have values T c > 4000 K, not low enough to prevent some 
emission from CELs, e specially in the [N II] lines ( see, for 
example, Figure 2 in ITsamis fc Pequignotil2005[ l. This 
emission will lower the value of T e ([N II])/T e ([0 III]) and 
hence can destroy the agreement between observations 
and models. There are many uncertainties related to 
the physical conditions of the hypothetical cold gas, and 
if its temperature is low enough, its emission will not 
contribute to CELs significantly. Hence, the observed 
values of T e ([N II])/T e ([0 III]) can be used as a further 
constraint on this kind of explanations. 

7. DISCUSSION 

Typical H II regions have a zone close to the ionizing 
star where ++ is the main ionization state of oxygen 
and a zone where + dominates (and where most of N + 
is also to be found). The stellar radiation field that de- 
termines this ionization structure heats the gas mainly 
through the ionization of neutral hydrogen. The excess 
energy carried away by the ejected electrons increases 
with distance from the ionizing star because the pho- 
toionization cross section of hydrogen is strongly peaked 
at the threshold energy of 13.6 eV, and the preferred 
absorption of photons with these energies hardens the 
ionizing radiation field. The gas cools mainly through 
recombination and through the collisional excitation of 
low-lying energy levels of abundant ions of heavy ele- 
ments, responsible for the production of CELs. At each 
point, an equilibrium temperature is reached through the 
balance of all the heating and cooling processes. In this 
way, the temperature and the metallicity of the gas are 
inextricably related. The hardening of the radiation field 
at increasing distances from the ionizing star and the 
very efficient cooling provided by the [O III] lines pro- 
duce a temperature gradient that translat es into the mea- 
sured values of T e ( [N II] ) /T e ( [O III]) > 1 (|Stasinskalll980t 
lOsterbrock fc Ferlandl [2006'). Finally, the total emission 
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in CELs like [O II] and [O III] lines is determined by both 
the ionization structure and the temperature structure of 
the nebula. 

The ADF found between the abundances implied by 
[O III] CELs and O II RLs not only introduces an un- 
certainty in the absolute value of the gas metallicity in 
observed nebulae; it might also require for its explana- 
tion significant modifications of the temperature struc- 
ture predicted by simple photoionization models. There- 
fore, a comparison between the temperature structures 
and metallicities measured in observed objects and pre- 
dicted by simple models can give indications of how much 
modification the models can accommodate. 

If we tried to build realistic models of the observed re- 
gions, we would need to specify several, very difficult to 
know, pieces of information such as the shape and lu- 
minosity of the ionizing radiation field emitted by each 
massive star in the region, the distance of all of these 
sources to the observed area, the amount and distri- 
bution of the intervening material between the ionizing 
sources and the observed area, and the density distri- 
bution of material in the observed area. Instead, our 
approach has been to construct grids of simple models 
and to see whether they reproduce the trends defined by 
complex real objects. To minimize the uncertainties re- 
lated to the measurement of line intensities, we use some 
of the best spectra available for H II regions. We com- 
pare the temperatures, T C ([N II]) and T e ([0 III]), and 
the oxygen abundances derived with the observed CELs 
with those derived following exactly the same procedure 
for grids of simple photoionization models with parame- 
ters that should encompass the ones that characterize the 
observed objects. The models are not tailored to repro- 
duce any object; they just use the best guesses in Cloudy 
for some parameters plus state-of-the-art model atmo- 
spheres. These models reproduce easily the observed 
values of the temperatures and the temperature struc- 
ture measured by the ratio T e ([N H])/T ([0 III]) when 
their input metallicities are close to those derived from 
CELs in the observed objects. 

We also consider some modifications to the simple 
models so that they also reproduce the observed O II 
RLs. We find that temperature fluctuations intro- 
duced by some extra heating mechanism will destroy 
the agreement found for the simple models by decreas- 
ing T e ([N II])/T e ([0 III]) to values below the observed 
ones. This problem might be fixed if the extra heating 
mechanism happens to be more efficient at larger dis- 
tances from the ionizing star. On the other hand, ex- 
planations of the ADF based on the addition of cold gas 
in metal-rich inclusions or ionized by X rays or cosmic 
rays, will only work if the temperature of this gas is so 
low that practically no CELs arise from the inclusions. 
None of these types of explanations of the ADF have been 
shown to work at the moment, and without more details 
on the mechanism behind the discrepancy, we can only 
provide these constraints that a successful explanation 
should meet. 

In view of the very good agreement found between the 
simple models and the observed objects, another expla- 
nation of the ADF in H II regions should be considered, 
namely, that the discrepancy is due to errors in the re- 
combination coefficients for O II. Some authors rule out 
this explanation because it would lead to similar values 
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of the ADF in all the objects, and though all H II re- 
gions and many planetary nebulae have ADF(0 ++ ) ~ 2, 
some plan etary nebulae have values that go up to 70 
(|Liul [2006). However, note that this does not rule out 
the possibility that the lowest values of the ADF, found 
in H II regions and many planetary nebulae, can be ex- 
plained in this way. Consider the Orion nebula, M42, 
the H II region with the best quality spectrum and 
hence the one where the measured ADF is more reli- 
able. Multiplet 1 implies ADF(0 ++ ) = 1.5 but other 
multiplets that should give reliable values of the ADF 
imply ADF{0 ++ ) = 1.04-1.77 (lEsteban et all 12001 . 
The same multiplets imply ADF(0 ++ ) = 1.5-2.3 in 
NGC 3576 (in the other nebulae only 1 to 3 reliable mul- 
tiplets were measured). If the recombination coefficients 
are underestimated by factors around or above 1.5 and if 
we also consider the uncertainties related to the measure- 
ment of these very weak lines, the derived ADFs could be 
easily explained in H II regions and those planetary neb- 
ulae with similar ADFs. A different explanation would 
only be needed for the handful of planetary nebulae that 
show high ADFs. 

On the other hand, the ADFs derived for other ions like 
C ++ , although much more uncertain than those derived 
for ++ , are usually larger than 1 and correlate roughly 
with ADF(0 ++ ) (see, e.g., Figure 18 in iWang fc Lull 
120071 for results on planetary nebulae). If the recombina- 
tion coefficients for ++ are underestimated, the same 
could happen with these ions. A revision of all these 
recombination coefficients would be valuable. 

Real nebulae are likely to have heating mechanisms 
besides those considered in models, and maybe some 
cold gas (due to shadows or ionized by alternative 
sources). This will introduce larger temperature fluc- 
tuations than those produced by simple models. In 
fact, the temperatures based on the Balmer or Paschen 
discontinuities of hydrogen seem to be significantly lower 
than those derived using CELs in s ome H II regions 
(see, e.g., iGarcia-Roias et all 12004 120061 ). However, 
in other H II regions, the temperatures implied by 
the hydroge n discontinuiti es and by CELs are fully 
consi stent (|Esteban et alj 120041 : IGarcfa-Roias et al] 
2005). The issue is far from being settled because the 
hydrogen temperatures are difficult to measure and 
have large uncertainties. But even if real objects have 
more temperature fluctuations than those predicted by 
models, most of the ADF could still arise from errors in 
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8. CONCLUSIONS 

We have constructed sets of simple photoionization 
models with parameters that encompass those of nine 
H II regions, chosen because their observed spectra are 
of high quality. Models and observations are analyzed in 
the same way to derive temperatures and oxygen abun- 
dances using the relative intensities of CELs and hydro- 
gen RLs. The good agreement found between observa- 
tions and models suggests that the models capture the 
main characteristics of the observed objects. In partic- 
ular, the agreement between the derived and predicted 
values of T e ([N II])/T e ([0 III]) implies that the temper- 
ature structure of the observed objects is reproduced by 
the models. Since the [N II] and [O III] diagnostic lines 
used to derive T e will react differently to changes in the 
temperature structure, this agreement places strong con- 
straints on the modifications or new ingredients that can 
be added to the models to explain the O II RLs. 

The most straightforward explanation of the ADF 
in H II regions (and those planetary nebulae with 
ADF(0 ++ ) ~ 2) requires errors by this amount in 
the recombination coefficients of O II RLs. However, 
we cannot rule out neither that there are relatively 
large temperature fluctuations due to some unknown 
heating mechanism whose efficiency increases with 
depth into the ionized cloud nor the presence of cold 
ionized gas if its temperature is low enough to suppress 
all emission in CELs. But whatever mechanism is 
invoked to explain the ADF, not only should it have a 
plausible origin, it should also satisfy the stringent con- 
straints imposed by the observed temperature structures. 
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